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Abstract 



Two methods for determining the potential (tp) around a discretely charged rod have been devised. 
The methods utilize the potential around the continuously charged rod (i/)) as the reference where i/; is 
determined by the Poisson-Boltzmann equation. The potential data are used to determine the theoretical 
radial distribution function (RDF) which is compared with MD simulation data. It is shown that the 
magnitude of the charge and size parameters very strongly affects the shape of the RDF's and consequently 
the thermodynamics. 

1 Introduction 

This work is motivated by the oscillations that occur at the polyelectrolyte (PE) radial distribution function 
(RDF) of a system simulation (see Fig. The PE chain in that figure are our model for the DNA polyelec- 
trolyte molecule. We use two kinds of PE models for the anionic bonded segment consisting of monomers of 
charge Zp. The first model has finite monomer charge Zp = — 12 where the distance between adjacent monomers 
is 6 = 20.4 A and the second model has z-p ~ —1 with b — 1.7 A. These models have the same line charge 
density Si = z-p/b. These models have the same radius and charge density as the B-DNA chain [H. The hard 
core radii of the DNA monomer is 8 A and the soft core radii is 2 A. For the temperatures set in the simulation, 
the effective DNA radius can be taken as 10 A because of the strength of the soft core repulsive potential which 
is effective over a distance of 2 A. Clearly we can observe oscillations at the RDF of the first model which does 
not appear in the second model. Thus the choice of the magnitude of the finite charge assigned to a model 
clearly affects the surrounding particle distribution. The degree of change of the RDF profiles with changing 
charge magnitude is very significant. In order to monitor these significant changes, here we focus on the simple 
model of an infinitely long charged rod where initially we focus charges being discrete and uniformly distributed 
along the rod. From the solution of the Poisson-Boltzmann equation (PBE) for the continuous (non-discrete) 
charge distribution of a charged rod, we can in at least two ways from the potential of the uniformly charged 
rod determine the potential field for a rod consisting of discrete charges. For a given (line) charge density Si, 
the charge magnitude and intermolecular distance are arbitrary. In order to achieve the above results, we first 
provide a theoretical analysis and continue with the numerical solution of the nonlinear Poisson-Boltzmann 
equation (PBE) for this rod model. To substantiate our theoretical analysis, we perform molecular dynamics 
(MD) simulation of this rod model of discrete charges which is surrounded by a salt solution. 
From the potential of a discretely charged rod, we determine theoretically the particle RDF about the rod. 

The method utilizes the Boltzmann factor (?f(r) = exp(— x^^-^) where ^° = g{ is the ratio of the number 

• • 

density of particle s at a distance r over the bulk number density; is the potential of average force [2]. 

This study is a first step in developing a model of calculating the general potential due to a discrete charge 
distribution enclosed by a boundary of arbitrary geometry. A principal motivation in attempting to provide a 
rigorous theory is due to the general tendency of workers to assume that the dimensions and parameters chosen 
from physical considerations for simulations is a fair representation of reality Q ; Our work shows that these 
choices are rather arbitrary and possibly inaccurate, which warrants a separate study of the connection between 
molecular size and the general potentials, since any one average potential of force ip(r) uniquely determines the 
RDF profiles. 



2 Theoretical Approach and Numerical Solution of PBE 

In the following we proceed to determine the non-uniform electric field distribution on a cylindrical Gaussian 
surface surrounding a charged rod. Then we relate it to the potential ip{r) and the particle distribution g{r). 



2.1 Nonuniform Electric Field Around a Charged Rod 



Fig. [T] shows an array of charges (closed circles) along the rod axis. The rod is extended to infinity to the 
left and right. A test charge c (open circle) lies at a perpendicular distance r to the rod axis. We start by 
determining the net electric field at point c due to the rod charges. For the moment, we assume no interference 
of counterion charges. We use bold letters to denote vectors and italic letters to represent its magnitude. In 
general, the first subscript in our variables refers to the left (value 1) half of the diagram and the subscript with 
value 2 refers to the right half of the same. The second subscript denotes the charge index along the rod in the 
appropriate direction indicated by the first subscript. 

The angle aj is the angle between the rod axis and a straight line from the point c to the left part of the 




Figure 1: The electric fields experienced at the point c (open circle) with a distance r normal from the polyelectrolyte 
rod axis. The closed circles represents the discrete charge distribution of the rod. 



rod charges. The index j starts from unity for the closest rod charge at the centre of the diagram and ends at 
infinity. Eij is the electric field at point c contributed by the j*'* charge on the left of the rod. Analogous to 
aj and Ei the angle /3j and electric field E2,j pertains to the angle and electric field for right part of the rod. 
In Fig. m if ri i denotes the distance between the test charge c and the nearest rod charge on the left side, p is 
the projection of the ri.i line to the rod x axis. Similarly, s is the projection of the line r2.i to the rod axis on 
the right hand side and s = {b — p). For later use, we define 

P,=P+{j-l)b (1) 
s,^ib-p) + {j~l)b (2) 

Below is the net electric field (in cgs unit) Ei due to the charges on the left. 

oo oo 

Ei-EEi.=Ep^r',,' (3) 

i=i j=i I'-J 



where r'lj is the unit vector from the j*^ charge to the test charge c and Zp is the magnitude of the discrete 
charges, or the total charge of a continuously charged rod over distance b. The prime coordinate are centred 
on test charge c, whereas the unprimed coordinates refer to the rod axis, where the zero for x is arbitrary, and 
the zero for the y coordinate is located on the rod. Components of r' , x' and y' , refer to orthogonal bases x' 
and y' that are parallel and perpendicular to the rod respectively. In mks units, we multiply by the coefficient 
l/Aire before the summation of the (|3]). We use cgs unit here for reasons of convenience. We can separate each 
Eij to their horizontal (x') and vertical (?/') components. Then Ei in terms of its x' and y' components are 
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Similarly for the electric field E2 due to the charges at the RHS of the rod 
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where the negative sign in the ^2,x expression is due to the reverse polarity of the x' axis. Since Ei and E2 act 
at the same point c, we can directly sum the net electric field for the x' and y' directions, which are 
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Then the magnitude of the resultant electric field acting at point c becomes 
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The ri j and r2.j expressions are (see Fig. [T]) 
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Fig. [5] depicts a continuously charged rod inside a closed cylindrical Gaussian surface. The figure uses {x',y') 
coordinate. Except for the electric field near the two ends of the rod, the direction of electric field is always 
perpendicular to the rod axis since the components parallel to the rod axis cancel everywhere. We define Zend 
as the length at the two ends of the rod where the electric field direction is not perpendicular to the rod axis. 
The Zend length is approximately finite depending on the rod charge and dimension. We define Loo as the length 
within the rod axis where the electric field directions are perpendicular to the rod axis. Since the rod has infinite 
length and Zend is finite, obviously Loo limits at infinity. From Gauss' law, the relationship between the electric 
field at the Gaussian surface and the total rod charge QRod is 0] 



E • dA = 47r Qr„ 



Breaking into the surface components of the closed cylindrical surface in Fig. [21 we obtain 
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Figure 2: Infinitely charged long rod (solid lines) inside a cylindrical closed surface (dashed lines). The continuous line 
charge density of the rod is 5i. 

where Ey is the electric field magnitude within the range Loo where the direction is perpendicular to the rod 
axis. i?y.cnd and i?x,ond are the electric field components from the Zend areas where the directions are parallel to 
the y' and x' axes respectively. 5i is the rod line charge density. Divide all terms in ([15]) by 27r r L^o , we obtain 
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Since the /end is finite, the numerators at the second and the third terms in the LHS of (fTB)l are also finite. Thus 

as Lqq ^ oo, 

(17) 



Ey = 



We define E = Ey as the electric field due to a continuously charged rod. Henceforth any barred symbols denote 
properties of a continuously charged rod model. 

If Zp is the total finite charge within an axial distance 5 of a continuously charged rod, then Zp = bSi and 
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If we define 
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then by dHHU) and p^ M^ . the following equalities hold 
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(22) 
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We define D as the electrostatic field distribution function. 

We utilize two kinds of discretely charged rod models shown in Fig. [3] with the same line and volume charge 
density. Figure |4] graphs the function Dy — Ey/ E for the two models. Fig. Hja for the rod model A whose finite 
charge —12 and Fig. Ulb for model B whose finite charge —1. Dy is determined numerically by extending j to a 
large number ('^ 5, 000) until no change is observed. From the figures, the electrostatic field distribution around 
the rod model B is rather uniform. For the rod model A the degree of variation of the curves for different r 
values is very obvious. The degree of variation diminishes as the perpendicular distance r increases. 

Figure [5] graphs the |I?2,(p)| = [iJjj/i?! for model A and model B. Like Dy, the variation of Dx appears only 
for model A and the variation vanishes as the perpendicular distance r increases. From Fig. |3]and[5l the values 
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Figure 3: Two charged rod models studied. Both model has the same line charged density zp /b 
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Figure 4: Function Dy{p) = Ey/E ([23} for: (Fig. a) rod model A {zp = -12 and b = 20.4A) and (Fig. b) rod model B 
(zp = — 1 and 6 = 1.7A). r is the perpendicular distance to the rod axis. 





Figure 5: Function |L'^(p)| = ^ for; (Fig. a) rod model A {zp = -12 and b = 20.4A) and (Fig. b) rod model 

B (zp = —1 and b = 1.7A). r is the perpendicular distance to the rod axis. 



of \Dx\ in general are much less than Dy. Since D = + D"^ (from (|19II21|) ) the contribution of Dy to D is 

very much larger compared to the contribution. 

In the following we propose a model when the counterion and salt ions are included in the system containing 
an infinitely charged rod. If the discrete rod charges along the rod axis have equal axial charge distance h, the 
ionic distribution about any two tangential plane that are equivalent by symmetry must be the same. We note 
that the D function i not affected by the counterion about the rod. Hence we may couple D to the effective 
dielectric constant that arise from the water and ionic distribution about the rod to relate E and i?,. 



2.2 Multigrid Method for Solving the PBE 



For some simple and symmetric shapes (e.g. plate, cylinder or sphere), the PBE expressions are well-known in 
their second order differential equation forms. The linearized PBE known as the Debye-Huckel approximation 
provides an analytical expressions that holds only for low potentials. As far as the second order differential 
form of the PBE is concerned, most solutions need to be found numerically. In this work we solve the PBE for 
both the continuous and discrete charged rod model. The potential for the continuous model will be used as a 
means to obtain the potential of the discretely charged rod by methods described in the two next subsections. 
For a continuous charged rod, the PBE is expressed as 

^ ^ s—1 

where 0(r) — eijj{r)lk-BT is the reduced potential and As = —e^ZsPs.ooH^k^T). The Neumann boundary 
conditions are obtained by applying Gauss' law at the rod-counterion contact distance r = = (op + rdon) 
and r = 7'c as follow 

d(j){r = ftp + 7'cion) ZpC^ 
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where Op is the rod radius and rcion is the counterion radius. Zp is the magnitude of the individual charge of 
the rod in the discrete case or the total charge contained within a distance b (for continuous rod). The cylinder 
cell radius Vc is determined by a non-local macroscopic rod charge density Pz^.^^^ (per unit volume), where 

-^Rod _ ^R°d /r,o^ 

Pznod - V ~ V ' ^ ' 

Simulation box Cylinder PBE box 

where ZrocI is the total charge of the rod in one simulation box; Vsimuiation box = ^Box"^ is the cubic simulation 
box volume with simulation box length Lbox which is set (see Section ^ . Since the length of the charged rod 
(with total charge Zroci) in one simulation box equals Lbox, the length of the cylindrical PBE per Z-nod charge 
also equals Lbox- Thus 

= r ^ 2 °^ (29) 
rc = ^J^Lbox (30) 

The second boundary condition (Eq. ([TT)) ) is applied for at least three reasons. First, at large distances from the 
rod axis, the potential is about constant implying zero potential gradient. Second, due to charge neutrality of 
the cylindrical box, theoretically the surface charge density equals zero at the box surface leading to zero electric 
field. Third, the systems studied in the theoretical PBE are substantiated by simulations utilizing the periodic 
boundary condition (PBC). If the same model like PBC is applied in PBE calculation, there is a mirroring 
shape of the potential between two neighboring boxes (in the PBE calculation) . We assume that our simulation 
box is large enough to warrant ((27|) . 

The multigrid method applied in this research uses the Full Approximation Storage (FAS) procedure described 
in [6] . The outline of the algorithm for our application is given in Appendix |A] based on a slight modification 
of the equations used in the code given in Q which was developed on the basis of reference p] . The way we 
differ follows. Reference uses another equation, and not Further, the code in Q is in 2-D with uniform 
and equal number of grid for both axes, whereas ours is both in 1 and 2-D where in 2-D the number of grid 
for both axes can be different; Q uses the Newton- Gauss- Seidel method for the smoothing procedures and we 
utilize the globally convergent Newton-Raphson and conjugate gradient method. Another difference is [Tl uses 
the Dirichlet boundary conditions, whereas we use the Neumann boundary conditions. Oberoi, et al., [8| and 
Hoist, et al., Q have also used multigrid methods to obtain the solution of the PBE. We had examined the 
uniqueness of the PBE solution by supplying different initial conditions that led to the same final solution. The 
convergence criteria is taken when the residual Euclidean norm of (I25[) is zero. Computationally, the iterations 
fluctuate about zero and we choose a tolerance number ('^ 2 x 10^^) of this norm that define our numerical 
solution. The physical variables used for our theoretical PBE calculations are exactly the same as for our MD 
simulations, so that the results from MD may be compared to the theoretical model. The detail of the physical 
variables of the MD systems and the theoretical models are given in Sectional 



2.3 Obtaining the Potential About a Discretely Charged Rod 



After examining the electric field around the rod axis, we relate it to the potential. Often, the scalar potential is 
the starting point for analysis and we adopt this route. We propose two methods to determine the non-uniform 
potential around a discretely charged rod. Both methods require the potential data around a continuously 
charged rod (0). The first method directly modifies the (j) by using the electric field distribution function D. 
The second method require us to perform 2D-PBE with boundary conditions specified in terms of E and D. 
For the first method (method I), we initially compute cf) of points (i, j) where index i refers to the x axis and j 
for the r axis, where the x grid distance are uniform, as with the r grid distances, but they need not be equal 
in length interval. The details of the solution for (f) was detailed in the previous section, where Eq. plays 
a pivotal rule. The solutions are determined for all (i,j) in (x, f) space. We then determine a relationship 
between (f) and (Eq. [38| . We then form a grid in (x,r) space which is exactly the same as (a;,r) space. 
We then compute (p over the grid using our solutions for <j), where for solution at grid (i,j) we require 4> 
at {(z,j-l)(z,j)(z,j + l)}. 

For method II, we also construct an grid over {x,r) space. We then determine the values of (j) a-t the 

grid boundaries. We also require the derivatives of 4> at this boundary which corresponds to the von-Neumann 
boundary conditions. We then use these boundary condition values to compute the 2-D PEE over a grid in 
(a;,r) space where Eq. ((42|l is solved. Subsection 15.11 details the consistency of these two approach. 



2.3.1 Method I: Modifying the Potential of a Continuously Charged Rod 

The starting point is the Poisson equation (mks unit) 

V.E(r') = ^, (31) 

where p{r') is the charge density at r' and e dielectric constant. The coordinate (r') where the electric field 
E(r') and charge density p{r') occur are outside the rod volume. We then assume the e at the medium outside 
the rod is constant (r' independent). If we denote p as the charge density around a discretely charged rod and 
similarly p for the continuous charged rod, then 
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(32) 



In cylindrical coordinate 



The y axis in Fig. [T] is identical with the coordinate r in cylindrical coordinate. We will use the r symbol to 
represent the coordinate aX 9 — Q 'm cylindrical coordinates. Therefore any arbitrary point r' pip is a function 
of r, 9 and x. The electric field and charge density for discrete rod model is (r, x) dependent and for the 
continuous model is (r) dependent. Then (|32p becomes 
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Substituting Er for Er in (P^HMj) . ((M)) becomes 
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In this work we approximate p{r') by the Boltzmann factor Ps(r') — ps^oo ZgC ^""^f'' ) where p^.oo and Zg are the 
bulk concentration and finite charge of species s. Since E = — V0 and E is only r dependent, then p7p becomes 
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We need the new potential (/)(a;, r) in the LHS. We propose to do numerical computation of ([38|) by replacing the 
derivatives with their numerical differencing form \^. For a second order numerical differencing, one possibility 
for IMl) is 



Dy{x„rj) + 

X ^ ] Ps, oo ZsB 



2h^ 



Dy{xi,rj + i)-Dy{xi,rj-i) _^ -Dx (aJj+i ,rj (a:^- 1 ,rj ) 



2hr 



2h^ 



4>{xi,rj+i)-24>{xi rj)+(t>{xi,rj^i) j <j)(a:i,rj.|-i)-0(a:i,rj.-i) 



-Zs4>[xi ,rj ) 



(39) 



For reasons that follow, only a 1-D root finding technique is required to solve (p9| . FiglHlb illustrates a sample 
point {xi,rj) and its four nearest neighboring points [(a;i+i, r^), (x^-i, r^), (x^, rj+i), (x;, rj_i)]. By supplying 
the and D values of the sample point {xi , rj ) and its neighboring point to the RHS of (15^ , the RHS of (15^ 
can be computed. We then determine the potential (j){xi,rj) for in the LHS by the method of root finding, 
i.e bisection or Newton-Raphson method. For the computation at the charged rod-solution interfaces, we can 
choose either the boundary conditions that obtain at the interface or reduce the order of numerical differencing 
into its first order. We use the boundary conditions specified in Subsection 12 . 21 for E values at boundaries. 
We close this subsection by giving an additional note on From Fig. |l]and[Sl when the perpendicular 
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Figure 6: (a). Electric field (E(r)) and potential {<j){r)) of a continuously charged rod as a function of coordinate r only. 

(b). An example of a central point {xi,rj) and its neighbor points for the <j){xi,rj) potential calculation with 
method I (Eq. [39]) 

distance r to the rod surface is far enough, the value of Dy is constant unity ( = 1) and the value of is 
constant zero ( = 0). It makes the term and in the square bracket of (pS)) equal zero. Then total value 
inside the square bracket of equals Dy = 1 leading to 0(a;, r) = 0(r). Thus at large distance, the potential 
around a discretely charged rod is equal to its continuous counterpart. 



2.3.2 Method II: Applying the 2-Dimensional PBE 

We propose another method to obtain the non-isotropic potential by including the electric field distribution 
function D in the numerical potential calculation of PBE. 
By dni) and we have 



Figure 7: ABCB plane. Left: Point A lies at a distance rp which is exactly perpendicular to a rod charge. Point B 
is perpendicular to the middle of two adjacent charges, point C and D lies at the cylinder box. Right: The 
plane ABCD is divided into grids for numerical computation. 
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The electric field over 9 coordinate is isotropic (see Fig. [7]), thus we can omit the terms. Then 
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Recall that [J| 
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Then by expressing the charge density p{r, x) by the Boltzmann factor, (|41l) becomes 
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where and Zi are the bulk density and finite charge of species i; s is the total species in the system except 
for the charged cylinder. 

Eq. (|42|) is a partial differential equation with independent variables r and x. In the x direction, the integration 
is performed over x G [0, 1/2&] because (j){r^x) is symmetric and periodic. We define a small rectangular plane 
ABCD (Fig El) with plane width AB = CD = 1/2&. The plane length BC = AD = r^- rp, where is the 
cylindrical box radius and rp is the closest contact distance between the rod axis and the counterion center. The 
corner point A is set to lie exactly perpendicular to the discrete rod charge with distance rp from the rod axis. 
Point B lies exactly on the perpendicular to the rod axis, pass through the mid-point of adjacent rod charges 
with distance rp from the rod axis. To solve the potential solution within the ABCD, we specify the Neumann 
boundary conditions at the edges of the ABCD plane. The potential derivatives which is the negative of the 
electric fields are specified by combining the electric field distribution function D and the electric fields of the 
continuous charged rod. Then the boundary conditions for (|42)) at the ABCD plane are 
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The electric field in the x direction of the boundary BC and AD equals zero because of electric field cancel- 
lation due to rod charge symmetry. E{r) is the electric field of a continuous charged rod model having the 
same dimension and line charge density as with the discretely charges rod. One possibility to obtain E{r) is to 
calculate the potential gradient of 0(r) by numerical differencing. Thus before we find the potential solution 
of a discretely charged rod, first we need to obtain the potential of the equivalent continuously charged rod to 
define the boundary conditions. 



3 Calculating the Mean Radial Distribution Function (RDF) from 
the Potential Data 

At this point we assume that we have been able to obtain the potential at any point around a charged rod. We 
then compare a potential-related property, which is the radial distribution function, obtained from simulation 
and calculation. Radial distribution function gs^si{f') = Ps,si{f') / Ps,oo measures ratio of the average to bulk 
density of particle s at radial distance r' from the central particle si. In this article, the central particle is the 
rod charge and the surrounding particles are the counterion and salt. As such we suppress the si subscript and 
only retain the s label. 

Since the potential from the cylindrical PBE is a function of the perpendicular distance r to the rod axis and 




Figure 8: Radial distance r' of a finite rod charge (closed circle). The rod axis is in the x direction. 

the lateral x distance from a fixed reference charge on the rod, we need a special treatment to obtain the mean 
RDF around a rod charge. 

Fig. [8]depicts a dashed-line hemisphere with radius r' and surface area A centered at a rod charge. The average 
number density at the hemisphere surface is 

P.ir')= ^sPsir.^^O)^A ^ ^^^^ 

where ps{r,x,0) is the number density at any point at the surface. The density is isotropic over 9 coordinate, 
thus we can exclude the 6 dependency. The surface area integration of (1471) follows the formula for a sphere. 
(|47)) then becomes 

_ 27rr' J^' p,ir,x)dx 

However, there is one difference between the calculation of bulk concentration of ions in PBE system and 
simulation. In PBE, the rod volume is not a space that is accessible to ions, thus the bulk density of ions is 
calculated by 

P^o^ = , (49) 

VBox — VRod 

where Ng is the total number of ion s, Vbox is the box volume and VRod is the rod volume. 

In simulation, the calculation of the ionic bulk density involves the whole box volume. So the bulk density of 
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where the definition of Ng and Vbox are similar with (|49p . 

The Ps(r') in (|48p is the source of the theoretical g{r') that will be compared to the simulation g{r'). Thus 

Ps{r') = pl%9{r') . (51) 

On the other hand, the local density ps{r,x) in (|48|) is the data from PBE calculation, so that 

Ps{r,x) = p™c^gf(r,a;) 

= exp(-z,0(r,x)) , (52) 

where gj has been defined for the PBE in the introduction. Substituting Ps{r') and ps{r,x) in (|511) and (|52|) to 
(l48t. we then have 



Pl%9{r') = ^ f p™- exp(-z,(/.(r,x))dx (53) 
Since p^"^ and pf™ are constant, ([55)1 reduces to 



„PBE 



£'(^') = - / exp(-z^(/)(r,x))da; . (54) 
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From and ([501) . since the total number of particle TVs is conserved in our scheme, we have p™^ (Vbox ~ 
VRod) = pf'^ Vbox- So dSll) becomes 

gir') ^ ^ f exp{-Zs(j){r,x))dx . (55) 

r Jo (Vbox-Vroci) 

If the rod volume is much smaller than the box volume (VRod ^ Vbox), the term Vbox/(Vbox ~ VRod) ~ 1 
Numerical integration of (j55p needs to be performed because there is no closed expression of 0(r, a;) which is 
numerically computed from (1551) and (|42p . We have used both in this presentation (i.e. Method I and Method 
II). The potential at any point {r,x) which is not exactly at one of the PBE solution points can be determined 
by interpolation. The interpolation method that is used in this work is either the polynomial or bicubic inter- 
polation 



4 Simulation Part 

We perform MD simulations for systems containing the rod model A and B (See Fig. [S]) in different salt concen- 
trations inside a cubic box. The rod dimension in the simulation box is made possible by making a cylindrical 
wall constraint. We fix the x and y of the rod axis in the middle of the simulation box (Xbox/S, iBox/2), where 
^Box is the box length. The rod sides are extended in the z direction and the periodic boundary condition 
applied will make both ends of the rod at infinity. We set the simulation box length to 408 A so that if the 
simulation is for rod model A, the number of discrete charges for the rod in the primary simulation cell (box) 
is 408/20.4 = 20. For the simulations of rod model B, the number of discrete charges for the rod per simulation 
box is 408/1.7 = 240. For both models, the total rod charges inside one simulation box is —240. Salt ions are 



modeled by spheres with radius 2 A where the charge is -1-1 for Na+ and —1 for Cl~ 13|,|lJ|. The counterions are 
chosen to have equal charge and radial dimensions as for Na+ from the salt. We use the Langevin thermostat 
to maintain the system temperature at 300 X for the whole MD run. The Bjcrrum length 7.13 is used as a 
parameter representing the continuum water solvent dielectric at this temperature. 

The purely repulsive Lennard- Jones potential is applied to control the short-range repulsive potential. We set 
the sum of the soft sphere radius (ctlj) for the LJ interaction between the rod surface and ion to 0.1 A implying 
that the sum of the hard sphere radii equals 11.9 A. The rod-ion soft sphere interaction radii is chosen to be 
very small to make the characteristic contact between ions and the rod surface as close as possible to the PBE 
boundary potential conditions. In PBE, the counterion cannot approach the rod axis at a distance that is 
smaller than the the sum of their effective radii (12 A). We use the ESPRcsSo package to run the MD sim- 
ulations . The RDF data are generated by the usual dumping procedure after the system attains equilibrium. 



5 Result and Discussion 



The purpose of this work is to examine the consequence of modelhng a bunch of charges as one single charge. 
This is frequently carried out in simulation studies that are interested in the physical forces that determine 
the particle distribution in complex systems |3|. In these studies, there is an underlying assumption that the 
simplification of molecular charge distribution in their simulations has a minimal, if not insignificant effect on 
the resulting general physical properties and distributions of the system. We explicitly show here that this tacit 
assumption needs to be reconsidered on the basis of Fig. [9] 

In order to clarify this situation, we present first our numerical methodology, and will use these techniques in 
subsequent writings where the PBE will be solved for a DNA polyelectrolyte with different monomer charge 
and size choices. Then, using our methodology, we can derive the RDF's theoretically and compare it with 
our simulation results. However, before this can be done, we need to ensure consistency of our method for 
some model systems, the most amenable being for the continuous and discretely charged rods. We therefore 
study in depth these systems in what follows. It may be added that this particular approach of the solution 
to this very important problem in biophysics concerning charge distribution profiles by using a continuously 
charged reference has not been attempted. The consistency of the result would imply that it is indeed a feasible 
methodology for all problems related to temperature-dependent charge density distributions. 
An example is when 12 (in number) independent charges in the DNA chain is simplified to a single charge, where 
the chain dimension is conserved. We had mentioned in the introduction concerning Fig. [9l where simulation of 
a linked spherical monomer representing the DNA chain results in the RDF being different when the monomer 
has a discrete charge —12 and —1 with the chain dimension and line charge density being conserved for both 
chains. Here, instead of working with a linked monomer to simulate an array of charge, we use the simpler 
model which is the infinitely charged rod for rod model A and B. The dimension of the rods are illustrated 
in Fig. [3] (for model A and B) where the rod line charge density are the same for the PEs whose RDFs are 
depicted in Fig. [HI 

In what follows we present the potential around the discretely charged rods (for model A and B), as a results 
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Figure 9: The RDF profile between the DNA polyelectrolyte charge and Na^ ions. The DNA monomer has charge —12 
(dotted line) and —1 (solid line). The details are discussed in the introduction. 

of the theoretical calculations utilizing method I and method II fsubsection l2.3.1l and l2.3.n) . From the potential 
data, we then construct the theoretical RDF by using the procedure given in Section [3] The theoretical RDFs 
are then compared to the RDFs from MD simulation. The theoretical potential, theoretical RDF, and simulation 
RDF are studied for both rod models A and B. 

5.1 Potential Data 

The equipotential contour lines around the charged rod from the calculations of method I and method II are 
provided (Fig. [TUl and [TT|) . From our numerical data, the potential profiles generated by method I and II are 
quantitatively the same within computational error suggesting that both methods are equivalent in the potential 
calculation. Method I is computational cheaper than method II because method I only performs a root-finding 
procedure at one particular point to obtain the potential at that point, whereas method II numerically calculates 
the 2-D PBE simultaneously for all the points within the rectangular grid. It may be added that the boundary 
conditions in method II gives some insight regarding the potential profile within the boundary whereas method 



I cannot provide such insight. 



5.1.1 Method I 





10 



10 



10 

x(A) 
(a) 



15 



20 



0.4 0.8 
x(A) 
(b) 



1.2 1.6 



Figure 10: Potential contour at a distance r perpendicular to a charged rod calculated by method I. Figure (a) for model 
A (zp = -12, 6 = 20.4 A) and Fig. (b) for model B {z-p = -1, 6 = 1.7 A) for NaCl sah concentration 10 mM. 



Figure [TOl shows the (j3{x,r) potential contour profile at a plane whose normal is perpendicular to the rod axis. 
The fig. [Tola for model A, and FigfTOlb for model B. The r axis is the perpendicular distance to the rod axis 
and the x axis is parallel to the rod axis. The range of x axis in Fig [10] equals the charge-charge distance b for 
each model. The two adjacent charges are at {x, r) — (0, 0) and {x, r) — {b, 0). The potential around rod model 
A varies and is symmetric. The largest absolute value of the potential around model A is at the nearest distance 
of approach of an ion to the discrete rod charge (at 12 A). In our simulation, the test charge has radius 2 A and 
the discrete charges of the rod all have radius 10 A leading to 12 A for the distance of closest approach. The 
potential differences disappear as the distance to the rod surface increases. For model B, there is no significant 
variation of the potential around the rod axis as x varies. 



5.1.2 Method II 
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Figure 11: Potential contour lines at a distance r perpendicular to a charged rod (model A) calculated by method II for 
different salt concentrations. 



Figure [TT] depicts the potential contours around the charged rod in different salt concentrations for rod model 
A (Fig. [3]) generated by method II calculations. The variables r and x are as for the previous subsection. 

Fig. [T2l depicts the difference between the potential calculated by method I and method II for the model A 
rod system in NaCI salt concentration 100 mM. It will be observed that the potential difference arc minute (of 




Figure 12: Potential difference (A<^ — (pn — 0j) from the potential calculation for system of rod model A using method 
I and II. Salt concentration 100 mM. 



the order of 1 % since (p w —3.6 at low r). We attribute the source of error to (i) numerical differencing and 
(ii) interpolation. 



5.2 Radial Distribution Functions 

From the theoretical potential data, the theoretical RDF is produced by the procedure of Section [3] Fig. [T51 
gives the RDF's of the Na+ ions about the rod charge of rod model A (Fig. [TSla) and model B (Fig. [T3lb) in 
different salt concentrations. 
In Fig. [T3]we observe some oscillations of the RDF for the simulation and theoretical RDF for rod model A 
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Figure 13: The theoretical rod charge— Na RDF (solid lines) and simulation RDF (dotted lines) for rod model A (Fig. 

A) and rod model B (Fig. B). The theoretical potentials about the discretely charged rod are calculated 
by method I (Subsection 12.3.1]) . The value 0.1; 10.0; 50.0 100.0 are the respective salt concentrations in 
millimolar. 



which are not observed for rod model B. We define oscillations have to be smooth and wavy slope of the graph. 
On the other hand the simulations also exhibits fluctuations in addition to the oscillations that occur in both. 
Little oscillations that occur at the simulation RDF of the rod model A has exactly the same period as the 
PBE RDF. The coincidence suggests that the principal cause of these oscillations in the simulation RDF is the 
nonuniform potential distribution parallel and perhaps perpendicular to the rod axis. For low salt concentra- 
tion, there is a disagreement between the height of the first peak of the simulation and calculated RDF, both 
for rod model A and B. This disagreement reduces as the salt concentration increases. The height of the first 
RDF peak between the simulation for model A and B is also different. This difference, which also appears for 
the RDF of DNA polyelectrolyte simulation in Fig. [HI indicates that the counterion distribution at the first 
layer about the charges of rod model A is denser than for the rod model B. 

The above discussions indicates that discrepancies can occur when simulating a bunch of charges as equivalent 
to one single charge. 



6 Conclusion 

We have shown that the potential profile for model A and B is quite different, implying that criteria based 
on theoretical modelling should be used in conjunction with MD simulations. Previously this aspect has been 
ignored in detail and elaborate studies that might open to question the validity of the results and conclusions. 
Method I and II generate identical results for the potential profile although the algorithms are entirely different. 
We can therefore exploit these methods differently where one method might be more tractable than the other 
for a particular systems. The theoretical and simulation RDF is in very good agreement, except at the first peak 
at low salt concentrations. We believe that the equations we solved are incomplete, especially in relation to the 
screening effects that are not taken into account. Finally it cannot be overemphasized that theoretical criteria 
must be adopted in choosing appropriate charge and size of the particles involved in the molecular system that 
is simulated. 
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A Multigrid method for solving the cyUndrical PBE 



To generate the potential data from (1251) . we apply the nested iteration multigrid method combined with the 
globally convergent Newton method (GCNM) for the smoothing process. In some cases, solving ([25)) with 
relaxation or shooting method cannot reach convergence. We also found that the GCNM offers better stability 
than the nonlinear conjugate gradient method for the smoothing process. The code is largely derived from Press, 
et al., 01, where we have modified the code, including replacing the original equation, the smoothing functions, 
and the boundary conditions (from Dirichlet to Neumann) in addition to expanding from 2-D calculation with 
equal number of grid for both axes to the 1-D and 2-D calculation where the number of grid at both axes can 
be different. 

Following the notation in we define C{u) (see Eq. (19.6.22) at Section 19 in [7]) as the source term for (I^St 
such that 

^ ^ i—l 

where u and x denote (j){r) and r respectively. 

The following is one of the modified portions of the file mgfas.c from (Press, et al., fi\). This C file is the driver 
for the nonlinear multigrid code. 



for(jj=j; jj>= 2; jj ) { /* downward stoke of the V */ 

relaxGCNRO; /* GCNR (Smoothing function) */ 

//relaxNLCG(); /* Non-linear Conjugate Gradient (Smoothing function)*/ 
} 

dtmpl= iu[l][2]; 

itmpl= renew_continuity(); /* solve the coarsest grid solution */ 

if(itmpl) printf(" Coarsest grid solved, mid point %f > %f', dtmpl, iu[l][2]); 

else iu[l][2] — dtmpl; /* no coarsest solution */ 
for(jj=2; jj i= j; jj-^-H) { /* Upward stroke of V */ 



relaxGCNRO; 

//relaxNLCGO; 

} 

update_neumannbc(); /* check and update at boundaries */ 

Figure [14] depicts the potential solutions of ([25]) in different salt concentrations using the above method. We 
notice that the potential rapidly converges to zero at higher NaCl concentrations, which is to be expected. The 
von Neumann boundary conditions are also satisfied in the vicinity of the boundary. 




Figure 14: The PBE potential solutions around a continuously charged rod calculated by multigrid method, r is the 
perpendicular distance to the rod axis. The rod model has continuous charge density 5 = Zc/b — 1/1. 7A 
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